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Abstract 

The flow of a viscous fluid along a curving pipe of fixed radius is 
driven by a pressure gradient. For a generally curving pipe it is the 
fluid flux which is constant along the pipe and so I correct fluid flow so- 
lutions of Dean (1928) and Topakoglu (1967) which assume constant 
pressure gradient. When the pipe is straight, the fluid adopts the 
parabolic velocity profile of Poiseuille flow; the spread of any contami- 
nant along the pipe is then described by the shear dispersion model of 
Taylor (1954) and its refinements by Mercer, Watt et al. (1994,1996). 
However, two conflicting effects occur in a generally curving pipe: 
viscosity skews the velocity profile which enhances the shear disper- 
sion; whereas in faster flow centrifugal effects establish secondary flows 
that reduce the shear dispersion. The two opposing effects cancel at 
a Reynolds number of about 15. Interestingly, the torsion of the pipe 
seems to have very little effect upon the flow or the dispersion, the cur- 
vature is by far the dominant influence. Lastly, curvature and torsion 
in the fluid flow significantly enhance the upstream tails of concentra- 
tion profiles in qualitative agreement with observations of dispersion 
in river flow. 
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Figure 1: perspective figure of a helical pipe as an example of the curv- 
ing pipes containing fluid flow and contaminant dispersion that is modelled 
herein. On the right is a schematic diagram of the orthogonal curvilinear 
coordinate system local to the curving centre line of the pipe: u points along 
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1 Introduction 

Consider the dispersion of a contaminant, with diffusivity n c , in the steady 
laminar flow, velocity field q, of a Newtonian fluid of density p and kinematic 
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viscosity v in an arbitrarily curving pipe of radius a such as the helical pipe 
shown in Figure 1. The flow is pumped by an overall pressure drop which 
maintains a fixed fluid flux along the pipe; that is, a constant mean velocity U 
is maintained irrespective of curvature. Non-dimensionalise quantities with 
respect to the pipe radius a, the cross-pipe diffusion time cl 2 /k c , and the 
reference pressure pvU/a. The Navier-Stokes and continuity equations for 
the steady, incompressible fluid flow then become 

Tlq ■ Vq = -Vp+ V 2 <? and V-q = 0, (1) 

where 1Z = all jv is the Reynolds number. The contaminant evolves accord- 
ing to the non-dimensional advection-diffusion equation 

^ + 7VVc = V 2 c, (2) 

where V = all / k c is the Peclet number. In typical liquids the Peclet number 
is much larger than the Reynolds number as their ratio, the Schmidt (or 
Prandtl) number a = V/1Z = v/k c , is normally large: approximately 10 3 for 
the diffusion of material in liquids [28, plll9, e.g.]; although only about 8 for 
the diffusion of heat [1, P597]; 1 whereas in typical gases the Schmidt number 
is roughly 1 and so the Peclet and Reynolds numbers are comparable. The 
analysis is interpreted with these two cases in mind of the Schmidt number 
being either 0(1000) or 0(1). 

Here we analyse the flow and dispersion in an arbitrarily curving cir- 
cular pipe. Most analysis of dispersion assumes a curved pipe is toroidal 
[28, 23, 14, 6, e.g.] and most experiments are performed in helical pipes 
[33, e.g.] (see further discussion by Berger et al. [3]). Neglecting molecular 
diffusivity, dispersion in toroidal flow has been characterised from analytic 
formulae by Ruthven [28] and numerical solutions by McConalogue [20] us- 
ing the residence time of different streamlines. Interestingly, using similar 
ideas Jones [15] deduced a regime of anomalous dispersion in a twisted but 
piecewise toroidal pipe. The fluid flow in helical pipes has been the sub- 
ject of recent analysis [12, 34, 19, 38, e.g.], whereas the flow in arbitrarily 
curving and twisting pipes has received little attention although Pedley [24] 
accounted for leading order effects of variable curvature but ignored torsion, 
and Gammack & Hydon [11] investigate flow in pipes with exponentially 
varying curvature and torsion. Here I take these analyses further by simulta- 
neously determining the fluid flow and the dispersion in arbitrarily curving 

lr Thc Prandtl number of water is 13.4 at 0°, 9.5 at 10°, 8.1 at 15°, 7.1 at 20° , 5.5 at 30° , 
4.3 at 40°, 3.0 at 60°, and 2.2 at 80°. Whereas the Prandtl number of air is 0.71 . 
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circular pipes. The analysis is restricted to parameter regimes where the 
fluid flow is laminar — because of the induced secondary circulation, laminar 
flow is stable to higher Reynolds numbers in a curving pipe [33, 3]. The 
results thus will also be important in the flow and dispersion in microfluidic 
channels [13, e.g.]. 

In §2.1 we establish an orthogonal coordinate system based upon the 
arbitrarily curving geometry of the circular pipe, although requiring that the 
curvature of the centre line varies smoothly. Let the centre line of the pipe be 
described by R(s) where s measures arclength along the centre line. Then a 
useful set of vectors for space in the vicinity of the pipe are the unit tangent 
u(s) = R' of the centre curve, unit normal p(s) and the unit binormal b(s) 
(see Figure 1). These vectors and the curvature k(s) and torsion r(s) of the 
pipe are all connected by the Frenet formulae [17, §8.7,e.g.] 



Throughout this article I use a dash to denote d/ds. In a thin pipe the non- 
dimensional velocity is approximately Poiseuille flow, u ~ 2(1— r 2 ). However, 
there are corrections of O(k) due to the curvature which are determined by 
solving the Navier-Stokes equations (1) for the fluid flow, see §2 where low 
order expressions agree with the careful analysis of the flow in helical pipes 
by Tuttle [34]. 

Centre manifold theory provides a rationale to form low-dimensional mod- 
els of dynamics as I elaborated in the overview [27]. Here we model the long- 
time evolution of the large scale dispersion of contaminant along the pipe. 
Models of "long-waves" or "slowly varying in space" dynamics are justified 
in the centre manifold approach by requiring resolved longitudinal spatial 
structures to have small wavenumber [25]. With this proviso, diffusion acts 
relatively quickly across the pipe to cause the contaminant concentration to 
be approximately constant in any cross-section: to leading order c ~ C(s,t) 
where C = c is an average over a pipe cross-section. The centre manifold 
analysis then systematically accounts for how variations along the pipe are 
affected by the varying velocity profile to disperse the contaminant. We thus 
deduce in §3, as a generalisation of Taylor's model [30, 31], the advection- 
diffusion model 



where for the case of a pipe of circular cross- section the effective diffusivity 



u' = up , 



p' = —ku + rb , b' = —rp . 



(3) 
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That is, the shear dispersion in a straight pipe, D = V 2 /48, is in a curved 
pipe modified by a factor approximately 



As found by others, secondary circulation caused by fluid inertia depresses 
the effective dispersion, by about k 2 [3.0(^/10) 2 + 0.38a 2 (ft/10) 4 ] , but only 
for Reynolds numbers sufficiently large. In slow viscous flow, pipe curvature 
actually enhances the effective dispersion by about 7.2k 2 — an effect that has 
apparently often been neglected [33, p317] in experimental determination of 
dispersion coefficients. 

• For dispersion in gas flow, with 0(1) Schmidt number a, the dispersion 
is depressed by secondary circulation only if TZ is greater than about 15 
as otherwise the viscous enhancement is stronger. Remarkably, if the 
Schmidt number a is small, less than about 0.5, inertial effects in the 
fluid flow enhance the effective dispersion for TZ greater than about 50. 

• For dispersion in liquids, with say (9(1000) Schmidt number a, the term 
in K 2 7?- 4 cr 2 dominates the other terms for Reynolds number greater than 
about 0.5. Hence, in liquids and due to secondary circulations due to 
inertia, I reaffirm the reduction in effective dispersion. 

Since the Dean number 2 D = this last is the shear dispersion cor- 

rection verified by Nunge [23] and Johnson [14] as being significant for D 2 cr 
greater than about 100. A limitation of the expression (5) is that it pre- 
dicts a physically unrealisable negative effective diffusivity for large enough 
Reynolds number. In most cases, Schmidt number a larger than 1, the term 
in K 2 a 2 TZ 2 dominates the correction. Hence to maintain a positive diffusion 
coefficient the Dean number D < 25 /^/a or equivalently crQ 2 < 650. The 
expression (5) is a low order approximation to the correct curve, describing 
the downwards curving shape on the left side of Figure 3 of Johnson [14], 
but needing the higher order corrections described in §3.4 to describe the 
dispersion coefficient at higher Dean number. 

For lower Reynolds number the qualitative deductions above vary from 
those of Nunge [23] because their dispersion coefficient is different, see their 

2 As discussed by Berger et al. [3, §2.1.1.2], there are various and conflicting definitions 
of the Dean number: Berger et al. recommended the use of D = 2^[kR. which I have 
adopted here. This Dean number could be viewed as t/kIZ for a Reynolds number based 
upon the pipe diameter rather than the radius that I have used. 
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equation (76). In particular, I predict that shear dispersion is frequently 
enhanced for gases, the reverse conclusion to that of Erdogan [9] and later 
Nunge [23, pp. 363, 375]. I argue that the differences occur because all previous 
work, based upon the fluid flow solutions of Dean [8, 7] and Topakoglu [32], 
have assumed that the pressure gradient is fixed in the expansion in cur- 
vature k — an adequate assumption for flow in a torus or helix where the 
curvature and the torsion are constant. But in a pipe of generally varying 
curvature and torsion, as developed here, it is the mean fluid flux which 
is fixed along the pipe, not the pressure gradient. 3 Since, for a constant 
pressure gradient the fluid flux varies with curvature and torsion — generally 
first decreasing with increasing torsion then later increasing with torsion, 
see Yamamoto [36] and its correction [37] — it follows that the mean pressure 
gradient (24) varies along a generally curving pipe. To check my computer 
algebra program (listed in Appendix A) I temporarily fixed the pressure gra- 
dient in a helical pipe and found the resulting dispersion coefficient to be 
exactly equivalent to that given by Nunge [23], equation (76), except that 
the one term in 1Z 2 <jh? (my k is their 1/A) is zero in my results — I conjecture 
theirs is in error in this term. Because of the requirement to fix the fluid flux 
I recommend the use of (5) instead of the earlier published models of shear 
dispersion. 

The error of 0(5) in the shear dispersion coefficient given by (5) encom- 
passes modifications due to torsion r and to variations in curvature n along 
the pipe: the parameter 5 corresponds to the parameter rj in Gammack & 
Hydon's analysis of exponentially varying pipes, k oc e vs . The torsion only 
affects the dispersion coefficient at 0(k 2 t 2 ), as see §3, and so does not ap- 
pear in (5). The effects of axial variations are reformulated as memory of the 
effective dispersion coefficient some distance upstream. Such memory effects 
in shear dispersion in varying channels were first recognised by Smith [29]. 

Using computer algebra it is also straightforward to determine both higher 
order corrections to the dispersion coefficient and high-order terms in the 
advection-diffusion equation itself. These terms may be either used to refine 
the approximations, or to give good estimates of the errors in a lower-order 
approximation. Earlier work by Mercer & Roberts [21] gave a sharp estimate 
for the limit of spatial resolution in a straight circular pipe. 



3 Evcn Gammack & Hydon [11, p363] appear to fix the pressure gradient in their expo- 
nentially varying pipes by requiring the second order pressure correction p2 oc sin £ , where 
£ is their angular variable, and so their pressure correction has zero mean. 
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2 The fluid flow 

The first task is to find the laminar viscous fluid flow in the curving pipe. 
The focus of the paper is the dispersion in the pipe by the flow, but there 
are enough interesting and relevant features in the fluid flow itself to be dis- 
cussed briefly here — in particular this section confirms aspects of my analysis 
by reproducing many results of other authors about steady laminar flow in 
curved and twisted pipes. 

We assume that the flow is steady as appropriate to flow driven by a con- 
stant pressure drop through a fixed pipe. However, to maintain everywhere 
constant fluid flux, the mean pressure gradient, p', varies with the curvature 
of the pipe as given in (24). 

2.1 The orthogonal curvilinear coordinate system 

Expressions for the flow are derived in an orthogonal curvilinear coordinate 
system matched to the geometry of the circular pipe. The orthogonal coor- 
dinate system has been used by Germano [12], Kao [16], Liu [19] and Ya- 
mamoto [36, 37] to investigate the structure of the fluid flow in helical pipes 
up to Dean numbers of 2, 000, and is well known in hydromagnetodynamics. 
One difference here is that we do not assume the pipe is helical, instead we 
allow arbitrary variations in the curvature and torsion of the pipe — the one 
important restriction is that the curvature and torsion must vary only slowly 
along the tube. Such slow variations along the pipe were also assumed by 
Murata [22] in their analysis of the flow in tubes bent sinusoidally in a plane, 
and Pedley [24] in a leading approximation to the effects of curvature. As 
shown schematically in Figure 1, positions in space are labeled by (s,r, •&) 
and have position vector 

r — R(s) + r cos dp + r sin# b, where 9 = •& + <p(s) (6) 

measures the angle from the plane of the normal p to the point r; thus 8 = 
corresponds to the inside of the local bend whereas 9 = ±ir corresponds to 
the outside. However, due to torsion in the shape of the pipe the reference 
plane of the orthogonal coordinate system must twist along the pipe by an 
amount <p(s) where 

d<P 
ds 
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The unit vectors and scale factors of this orthogonal coordinate system are 
then 

h s = 1 — kt cos 9 , e s = u , 

h r — 1 , e r = cos # p + sin # b , (7) 

h$ = r , e$ = — sin # p + cos 9 b . 

Note that all expressions for fluid and concentration fields are written in 
terms of 9, the angle relative to the local direction of curvature of the bent 
pipe — because it is this angle that primarily determines the shape of the local 
fields — but all equations are written in terms of the angular coordinate in the 
orthogonal system, namely remember that 9 varies with $ and s according 
to (6). Observe the scale factors are all positive provided < r < 1/k 
and so the coordinate system is well defined for unit radius pipes provided 
the non-dimensional centre line curvature k < 1 . Let the velocity field, 
with components the axial velocity u, the radial velocity v , and the angular 
velocity w, be denoted by 

q = ue s + ve r + we$ . 

Then, noting it is convenient to compute the viscous dissipation term via the 
vorticity (as does Tuttle [34, p548]), 

V 2 q = —V x uj , u> = V x q , 

standard formulae apply for computing components of the Navier-Stokes 
equations (1) [1, Appendix B, e.g.]: 

1 / d(rw) dv\ 
r \ dr 8$ ) 



rh s \ d"d ds 
1 fdv d(h 8 u) \ 



h s \ds dr 
q 1 dp 1 ( d{ruj$) du T 



h s ds r \ dr dfi 

( u du du w du uv dh s uw dh s , 

\h s ds dr r dfi h s dr rh s dd 

dp 1 f d(h s uj s ) d{ruji t 



dr rh s \ dfi ds 

+ 7Z ( U ^ V + ^ V + W ^ V w2 ^ ^ s ) 
V h* ds dr r d-d r h* dr 
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h s ds dr r dd rh s dd r 



u dw dw w dw u 2 dh s vw 




(13) 







1 / d(ru) d(rh s v) d(h s w)\ 
rh s \ ds dr dfi J 



(14) 



These are solved with a fixed fluid flux and with zero velocity on the pipe 
walls: u = v = w = 0onr = l. The computer algebra program in Ap- 
pendix A solves these equations iteratively. 

There are some subtleties in the geometry of the coordinate system. As 
discussed by Zabielski [38, §2.2], observe that because of the twist in a helical 
pipe the axial unit vector u is not everywhere tangent to the lines of helical 
symmetry — the s-coordinate curves are not curves of helical symmetry. Thus 
be careful in interpreting cross-flow velocities v and w because in one view 
they will involve a small component of the relatively large velocity along the 
lines of helical symmetry. In an alternative presented by Tuttle [34], the 
twist in the coordinate system caused by torsion generates an effect similar 
to that caused by a coordinate system rotating in time. However, here we 
consider flow in a generally curving pipe with no large scale symmetry, so 
the only definite longitudinal direction is the local unit vector u and we thus 
discuss v and w as cross-flow velocities, as does Gammack & Hydon [11]. 
Similarly, in helical symmetry one cannot find a cross-section plane normal 
to the lines of helical symmetry [38, p300] so an arbitrary decision is needed. 
As is conventional for helical pipes and as simplest for generally curving 
pipes, we conventionally take a cross-section to be normal to the centreline 
of the pipe. In these cross-sections the pipe is circular. 



Solving the fluid equations using the computer algebra program in Ap- 
pendix A I deduce the Stokes flow field, 71 = 0, is 



2.2 Slow Stokes flow 



w 



u 



V 



(l - r 2 ) [2 + k\t cos 9 + K 2 |r 2 cos 29 - k 2 ||(1 - 3r 2 )] 

+ 0(k 3 ,5 2 ,TZ), 
| ( 1 — r 2 ) 2 [cos 9 k + sin 9 kt\ 

+ (l -r 2 ) 2 [(38 + 43 cos 29)kk +43 sin 29 k 2 t] 
+ 0(k 3 ,5 2 ,TZ), 

|(l - r 2 ) (2 - r 2 ) [- sin 9 k' + cos 9 kt] 



(15) 



(16) 
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Figure 2: (a) contours of axial velocity u of the viscously dominated Stokes 
flow in a helical pipe with curvature k = 0.8 (chosen so large to accentuate the 
modifications); (b) corresponding torsion induced cross-pipe fluid velocities 
to leading order in the torsion. The plots are evaluated from the asymptotic 
solution with errors 0(k 5 ). 



+ ±r (1 -r 2 ) [(43- 29r 2 )(cos2#K 2 r-sin2#KK') + (6 - 2r 2 )n 2 r] 
+ 0(k\5\TZ), (17) 

p' = -8 + \k 2 + 0(^,5 2 ,n), (18) 

p = p — |r (l — 3r 2 ) [cos 6* k + sin 6 kt] 

- JL (5 + 4 r 2 _ 21r 4 ) kk' - ±r 2 (9 - 26r 2 ) [cos 29 kk' + sin 29 k 2 t] 
+ 0(k 3 ,5 2 ,1Z), (19) 

where 5 is used to denote the order of magnitude of derivatives of the quan- 
tities varying slowly along the pipe. For example, k' and r = —(f)' are 
thus 0{5). 

See that, for example, the Stokes flow in a torus (k = const and r = 0) is 
simply one of axial flow, see Figure 2(a), in an adjusted mean pressure gradi- 
ent as all other components vanish. The axial velocity maximum is shifted to 
the inside of the curve (to the right in Figure 2(a)) and is increased slightly. 
In contrast to flows at significant Reynolds number, the pressure gradient 
around a curve is less than that in a straight pipe presumably because the 
bulk of the fluid travels a shorter path than the centre line — this agrees with 
Larrain [18] who used computer algebra to also find high order approxima- 
tions to the flow in a coiled pipe. 

The cross-pipe velocities in a helical pipe are indicated in Figure 2(b) 
where the torsion induces velocities proportional to those shown in the fig- 
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(a) 



(b) 




Figure 3: (a) axial velocity u contours of the inertial corrections to Poiseuille 
flow in a helical pipe with curvature k = 0.5 to leading order in the torsion r 
of the helix; (b) corresponding inertia induced cross-pipe fluid velocities inde- 
pendent of the torsion. The plots are evaluated from an asymptotic solution 
with errors 0(k 5 ). 

ure; the generally upwards velocity matching the upwards twist of positive 
torsion. Observe that torsion, r, and variations in curvature, k', only affect 
the cross-stream velocities and do not influence the axial velocity u to this 
order. Conversely, observe that this Stokes flow does not have cross-pipe 
circulation — the strong viscosity eliminates inertia. Instead the curvature of 
the pipe just skews and alters the velocity field. For curvature k ^ the 
maximum of the axial velocity u increases and moves towards the inner wall 
of the pipe. This last effect, though seemingly small even for the large cur- 
vature of k = 0.8 used in Figure 2, is enough to have a strong influence on 
the shear dispersion as seen in (5). 

2.3 Laminar flow at finite Reynolds number 

Incorporating the terms representing the advection of fluid momentum into 
the computer algebra program of Appendix A leads to effects parametrised 
by the Reynolds number 1Z. I find that the fluid fields given previously 
in (15-19) are modified by the addition of the following terms: 



u 



••• + 



144 



r 



(l - r 2 ) (29 + 5r 2 - 3r 4 ) [cos 9 k' + sin 9 act] 



n 2 



2 



1440 



r 



(l - r 2 ) (19 - 21r 2 + 9r 4 - r 6 ) cos 9 k 
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7?s 

+ r (l - r 2 ) (2969 - 4381r 2 + 3249r 4 - 1301r 6 

1814400 v ; v 



+ 274r 8 - 20r 10 ) [cos 9 k! + sin 9 kt] + 0{k 2 , 5 2 ) , (20) 
K 
72 



- ^ (l - r 2 ) 2 (4 - r 2 ) cos « 



7e 2 

+ (1 - r 2 ) " (13 - 15r 2 + 7r 4 - r 6 ) [cos k' + sin 9 kt] 

8640 V ; V n J 

+ 0(k 2 ,5 2 ), (21) 
■ ■ ■ + ^ (1 - r 2 ) (4 - 23r 2 + 7r 4 ) sin « 

(1 - r 2 ) (13 - 224r 2 + 266r 4 - 124r 6 



w 



8640 

+ 17r 8 ) [- sin 9 k' + cos 9 kt] + 0(k 2 , 5 2 ) , (22) 

p = ■ ■ ■ r (9 — 6r 2 + 2r 4 ) cos 9 k 

3 

+ (101 - 120r 2 + 90r 4 - 30r 6 + 3r 8 ) [cos 9k + sin kt] 

+ 0(k 2 ,5 2 ), (23) 

where "• • •" denote the terms already given for Stokes flow in (15-19). The 
modifications to the cross-pipe velocities that are proportional to 1Zk{= 
® 2 /41l), plotted in Fi gure 3, agree with those of the Dean flow as used by 
Johnson [14, p330], in their work on dispersion. The cross-pipe velocity field 
exhibits circulation across the pipe induced by the pipe curvature because of 
fluid inertia. The term in the axial velocity u proportional to kR 2 = £) 2 /4 
also agrees with that of Johnson & Kamm. These terms in the veloc- 
ity fields are those previously found for the "loosely coiled limit" [3, p467] 
when curvature k is negligible by itself but the Dean number D = 2tJk1Z 
is significant. The above expressions appear to agree precisely with the ex- 
pressions (58-61) carefully obtained by Tuttle [34] for low Reynolds number 
flow in a helical pipe — the only differences lie in various factors of two due to 
the different non-dimensionalisation and because my angular velocity w is in 
a "space-centred" coordinate system whereas Tuttle's $ is "body centred". 4 

4 Although the components of the velocity proportional to curvature k and the leading 
order terms in TZk' reduce to those of Murata [22, Eqn. (21-22)] , in their case of a sinusoidal 
centreline, the terms in 1Z 2 are different in detail to those of Murata [22, Eqn. (22)], as 
is the pressure. Terms in curvature n in the above velocity field agree with those of 
Pedley [24, Eqn. (4.13)], and terms in the gradient k' with highest power of Reynolds 
number 1Z also agree [24, (4.18-19)] except for the axial velocity u. The above expressions 
also agree with the small Dean number expansion derived by Kao [16, p341], for flow in a 
helix except for his wi which does not match my expression (20) for u. 



Tony Roberts, February 8, 2008 



2 The fluid flow 



13 




.4 1 1 . 1 1 

2 4 6 8 10 

Re 

Figure 4: coefficients in the expression (24) for the mean pressure gradient 
as a function of the centreline curvature k and its gradient k' . 

Lastly, the components of the physical fields (15-23) involving sin# and cos# 
cross-sectional structures agree with those of Hammack & Hydon [11, p363], 
but their formulae have none of the components in sin 29, cos 29, nor a mod- 
ification to the mean pressure gradient as in (24). Observe in all of the 
formulae (16-23) that torsion coupled to curvature, kt, appears to have the 
same effect as longitudinal gradients of curvature, k' , but the fluid fields 
are rotated in angle by 90°. The above expressions are the first to combine 
torsion and general variations in curvature. 

With fixed fluid flux, the mean pressure gradient, to one higher order in 
curvature than the fields above, is 

rl = a 1 (1 11 7? 2 1541 y 4 ) Jr 2 

r ~ V 6 540 /V 32659200 /V I 

+ (m + iHo^ 3 + iiSoo^ 5 ) + * 2 ) • (24) 

See from the coefficients plotted in Figure 4, that although the pressure gra- 
dient is lessened by curvature for low Reynolds number, for Reynolds number 
approximately 1Z > 3 there is an increased pressure gradient loss in a curving 
pipe. This is attributed to the greater mixing caused by the induced cross- 
pipe circulation. Also see that a region of tightening curvature, increasing k 2 , 
has a lesser drop in pressure gradient relative to that for a toroidal pipe, 
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whereas conversely a lessening in the curvature, decreasing k 2 , has a higher 
drop in the pressure gradient. This effect is interpreted as a "memory" in 
the mean pressure gradient of the upstream conditions: retaining just the 
highest order terms in 1Z the mean pressure gradient 

ft -8 - 1541 7? 4 ;r 2 4- 6191 7? 5 2*V 

r ° 32659200 /V ' 821145600 /V 

D 4 \ (2F>) 



1541 



32659200x16 I 5-4337^/271216 ' 

where the evaluation of the Dean number at an effective distance upstream 
of approximately 71/ 6 pipe radii seems to show the typical distance necessary 
for the fluid flow to develop in order to accord with the curvature of the pipe. 
This agrees qualitatively with experiments on a pipe with a finite bend as 
discussed by Berger et al. [3, p494], where the influence of the bend on the 
mean pressure gradient extends far downstream. This upstream memory also 
matches nicely with the commonly quoted distance, l s ~ \oR- [10] and [3, 
p488], required for flow entering a straight pipe to become fully developed, 
and with the observation by Murata [22, §4], that the flow in a sinusoidally 
bent pipe has a lag in its adaptation to the local conditions, the lag increasing 
with increasing Reynolds number. 

The above formulae for both viscous and inertia effects of curvature and 
torsion are only of low-order. Computer algebra computes the velocity field to 
as high an order as is necessary for the demands of modelling the dispersion 
in the pipe, described in the next section. Van Dyke [3, p475] has shown 
the asymptotic expansions of the fluid flow field converge for Dean number 
D < 96.8/4a/2 = 17.1 (for negligible k but finite 3)). However, I have not 
explored this issue as here we are primarily concerned with the advection- 
diffusion model (2) of the dispersion and its asymptotic approximations. 



3 Advection-dispersion along the curving 
pipe 

Having determined the fluid flow within the pipe, I now address the advection 
and longitudinal dispersion within the pipe. We solve the advect ion-diffusion 
equation (2) for the evolving concentration c(s, r, t) of contaminant within 
the fluid. Writing the concentration in terms of the cross-pipe average C(s, t) 
and its derivatives I use centre manifold techniques to construct the Taylor 
model (4), and its higher order generalisations, of the advection-dispersion 
in the bent and twisted pipe. 
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3.1 The concentration field within the pipe 



Using the coordinate system described in §2.1, the advection-diffusion equa- 
tion (2) for the evolution of the concentration of the contaminant is 



dc ( u dc dc w dc 
dt V ds dr r dfi 



1 



d ( r dc 



d 



dc 



ds \ h s ds ) ~^ dr \ s dr / d$ \ r d$ 



d I ' h q dc 



(26) 



This is solved with no flux through the circular walls of the pipe: dc/dr = 
on r = 1 . 

The computer algebra program (Appendix A) simultaneously determines 
from the contaminant conservation equation (26) the dynamics on the low- 
dimensional centre manifold, namely the Taylor model (4). Over a cross-pipe 
diffusion time the concentration field evolves to be, for example, approxi- 
mately 



6 V ds24^ 2 

n 2 



6r 2 + 3r 4 ) +V^kcos6 
' ds 



- (4r - 3r 3 + r 5 ) 
6 



172800 
all 2 



256r - 285H + 200r 5 - 75r 7 + 15r b 



,ii 



) 



H (68r - 120r 3 + 130r 5 - 75r 7 + 21r 9 - 2r n ) 

34560 v ; 

0(k 2 ,5 2 ). 



(27) 



To this order neither torsion nor gradients of curvature affect the concen- 
tration field within the pipe, but this is not surprising as one order of 5 is 
counted in dC / ds leaving no scope for derivatives of the curvature k to be in- 
volved in the above terms. Expressions for the concentration field to the next 
order in either the curvature k or longitudinal derivatives 5 are algebraically 
formidable and are not recorded. 

From such expressions, the computer algebra determines the mean flux 
of contaminant through a cross-section of the pipe: 



F(C, 



uc — 



]_dc 

h q ds 



where the overbar denotes the average over a cross-section. Then by con- 
servation of contaminant, the model for the evolution of the contaminant 
is known to follow C t = —F s . Using the flux F like this I determine the 
right-hand side of the model to one order higher in axial derivatives, 5, than 
would otherwise be possible. 
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3.2 Arbitrary curvature causes upstream memory 

The error of 0{8) in the shear dispersion coefficient given by (5) encompasses 
modifications due to both torsion r and to variations along the pipe in the 
curvature k. The torsion only affects the dispersion coefficient at 0(k 2 t 2 ) 
and so does not show up in (5). 

Variations in curvature along the pipe (k' ^ 0) cause the effective dis- 
persion coefficient to become, using 91 to denote the scaled Reynolds num- 
ber 7^-/10 and remembering V = TZa, 

D = 1 + -k 2 + — (kk'V - — ( k 2 t 2 + k' 2 ) (28a) 
4 48 V 1 96 V / v ' 

V 2 

, ' |"i , Jl ( 64225 o>4 2 , 2995 $vv4 36335 q}2 , 863 \1 fOSM 
48 ^ (~ 171072 ^ ^ + 24192^ ~ 12096^ + 12o)J ( 28b ) 

, T> KK ' IY»6 I 9050586625 4 , 246093875 3 _ 6234774125 2 _ 1760495125 \ 
L V 26900729856 " r 1630347264 53801459712 40351094784 / 

_|_ m 4 f 1068925 3 , 33738035 2 , 2310385 \ 

-T- ~n \ 2322432 ° *~ 20901888 4478976 / 

_|_ m 2 f 383695 2 , 19465 , 985 \ 

" l ~ V 96768 12096 12096/ 

- §] (28c) 

I -p2/ .2 2 _,_ /2n r~6 / 2542365125 4 _ 1039029345155 2 _ 5542735225 \ 
-VI ^/t f /t _/ |y\ l,5 86 9 2 50i504 O 126541033242624° 3515028701184/ 



, <v>4 / 14791164485 „1 , 739414405 



+ <n 2 ( 



33108590592 ~ 1 11036196864/ 
7181 2 , 7619671 \ 



9072^ 1 41803776/ 

- I§] (28d) 

-\-VW(kk'Y T^R 6 f 2219783253125 5 _ 3007270625375 4 , 7670650920025 3 
■VI J\\hl% j yj\ y 3012881743872 9038645231616 63270516621312 

, 2879131496575 2 , 396673566125 \ 
*~ 23726443732992 15817629155328 ) 

, o>4 / 318184675 4 _ 3159503125 3 _ 9765145925 2 _ 2663242675 \ 
^ V 306561024° " 752467968 16554295296 5518098432°/ 

2 / 40508065 ^3 3081595 Jl 3844625 ^ 16835 A 



_i_ <v>2 1 4U5U8U65 6 _ 3081595 2 _ 3844625 _ 



4644864 " 1548288 2612736 373248/ 



+ fH- - m + 0(^,S 3 ) , (28e) 
In this large but comprehensive expression observe: 

(28a) gives the molecular diffusivity along the pipe in the presence of the 
bending and twisting of the pipe when there is no flow; 

(28b) gives the usual shear enhanced dispersion in a straight pipe, P 2 /48, 
modified by the leading order (quadratic) effects of pipe curvature 
(these were the terms of the shear dispersion discussed in the Intro- 
duction (5); 

(28c) gives the leading order effects on the dispersion due to variations in 
curvature along the pipe; 
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(28d) gives the leading order effect of torsion on the dispersion, namely 
quadratic but moderated by the multiplication by k 2 ; 

(28d— 28e) through k' 2 and k" terms, gives the second order effects of the 
variations in curvature. 

It is intriguing to see that the effects of torsion and second order gradi- 
ents of curvature factorise as shown in (28d-28e). I suggest the reason 
for this factorisation is due to two effects: firstly, upstream "memory" 
of the dispersion, to be discussed later, involves 




which may explain the appearance of the combination (kk')'] and sec- 
ondly, curvature gradients and torsion, k' and kt respectively, both 
create the same but orthogonal structures in the fluid flow as com- 
mented after (20-23). 

For large Schmidt number er (typical for material dispersion in liquids) 
there are two distinguished limits of the above expression for the effective 
dispersion coefficient, the second being a subset of the first. 

• Firstly, for large Schmidt number a the highest powers of a domi- 
nate. However, in various subexpressions they appear in combina- 
tion with the Reynolds number TZ = 10D\. Thus there is a distin- 
guished limit with large a and small 1Z in which 1Z 2 <j is of order 1 . 
In terms of the magnitude S of the slow axial variations, an appropri- 
ate scaling is that the fluid flow is slow, 1Z ~ 8, the Schmidt num- 
ber large enough, a ~ l/<5 2 , so that the Peclet number is also large, 
V ~ 1/5, then the effective diffusion coefficient is large, D ~ 1/5 2 . 
Using these orders of magnitude, introducing the order 1 parameter 
a = TZV /100 = TZ 2 a/100 and evaluating fractions, the leading order 
terms in the dispersion coefficient (28) are 

D w .02083 P 2 + (.1498- .007821 a 2 )V 2 n 2 

+ (-.03965 - .004603 a + .003364 a 2 )PW 

+ (-.007916 + .0004332 a 2 )V 4 (K 2 r 2 + k' 2 ) 

+ (.008721 + .001038 a - .0007368 a 2 )V 4 (Kn')' . (29) 

• Secondly, for a typical Schmidt number a bigger than 10 3 or so, and 
for any flow with Reynolds number 1Z bigger than about 1, then the 
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parameter a will be bigger than about 10 and the above expression (29) 
will be dominated by the quadratic powers in a. That is, the dispersion 
coefficient 



We noted in (25) that the mean pressure gradient in the fluid flow at any 
location was appropriate to the curvature some distance upstream. Similar 
memory effects are seen in the dispersion coefficient. The subexpression 
—0.7821 k 2 + 0.3364 Vkk! appearing in the first line of (30) is equivalent to 
simply —0.7821 k 2 evaluated at a distance £ = 0.2151 V upstream from any 
particular location. I do not attempt to complicate this memory effect any 
further by trying to include the second order term ( there 
plethora of possibilities, but for the purposes of discussion I assume both kk' 
and (kk')' terms are attributable to upstream memory. The ratio of the 
coefficients of kk' and k 2 in (29) similarly quantify the upstream memory for 
low Reynolds number flows as shown in Figure 5. Such memory effects in 
shear dispersion in varying channels were first recognised by Smith [29]. 

From the coefficient approximations (29) and (30) see that torsion and 
curvature gradients generally enhance dispersion along the pipe except for 
low Reynolds numbers, TZ 2 a < 427. 4 , when they make the dispersion coeffi- 
cient smaller. However, the effect torsion has upon the dispersion coefficient 
seems small because not only is the effect quadratic in the torsion r, it is 
also ameliorated by the multiplication by the curvature squared. However, 
ignoring the kk' terms and noting that k! = /■c(logK)', I write (30) as 



This suggests that torsion, or proportional gradients of curvature, greater 
than about 4/V may cause the dispersion coefficient D to increase with 
curvature k, instead of decreasing. That torsion could eliminate the increased 
mixing due to secondary circulations seems unlikely so I predict higher order 
terms in the torsion r would limit its influence on the dispersion. 




(30) 




(31) 



Tony Roberts, February 8, 2008 



3 Advection-dispersion along the curving pipe 19 

0.3 
0.28 
0.26 
0.24 

CD 

S>0.22 
0.2 
0.18 
0.16 

50 100 150 200 250 300 350 400 

Re 2 Sc 

Figure 5: in a generally curving pipe the effective dispersion has the value 
appropriate to the curvature a distance £ upstream. 

For the dispersion of heat in water with a Prandtl number of er = 8.1 
at 15°C the dispersion coefficient (28) reduces to 

D = 1.3669 ft 2 + 1 

+ k 2 (-.003350 ft 6 - .04106 ft 4 + 9.830 ft 2 + .25) 

+ kk' (.01232ft 7 - .1090 ft 5 - 2.01 ft 3 - 3.291 ft) 

+ («V + k' 2 ) (.01220 ft 8 + .1928 ft 6 - 33.95 ft 4 - 6.356 ft 2 - .07292) 

+ {kk') 1 (-.02191ft 8 + .1777ft 6 + 36.39 ft 4 + 7.602 ft 2 + .1458) . (32) 

For Reynolds number ft > 7 the highest powers in ft dominate; see Figure 6 
for the dependence on smaller ft. Observe: for Reynolds number ft > 6.95 
curvature enhances the dispersion of heat and vice-versa; whereas for ft > 
6.74 torsion and curvature gradients reduce the dispersion and vice- versa. 

For the dispersion of heat in air with a Prandtl number of a = 0.71 
at 15°C and recalling that 9^ = ft/10, the dispersion coefficient (28) reduces 
to 

D= 1.080*R 2 + 1 
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Figure 6: coefficients of the dispersion coefficient D for water at 15°C given 
by (32). For higher Reynolds number the highest powers in TZ dominate the 
coefficients. 

+ k 2 (-.07649 <R 6 - 3.244 m 4 + 7.767 <K 2 + .25) 

+ kk! (.3979 9t 7 + 7.462 ^ 5 - 5.871 9\ 3 - 2.925 9t) 

+ (kV + k' 2 ) (.3011 + 15.48 - 11.82 5H 4 - 5.022 ^H 2 - .07292) 

+ {kk')' (-.7615 - 13.98 9* 6 + 8.055 + 5.282 + .1458) . (33) 

For Reynolds number 1Z > 15 the highest powers in 1Z dominate; see Figure 7 
for the dependence on smaller TZ. Observe: for Reynolds number 7Z > 15.18 
curvature enhances the dispersion of heat and vice-versa; whereas for TZ > 
10.35 torsion and curvature gradients reduce the dispersion and vice-versa. 

3.3 Skewness is very sensitive to curvature 

Computer algebra straightforwardly determines high order terms in the advec- 
t ion- diffusion equation (2). Chatwin [4] investigated the relatively slow ap- 
proach to normality in shear dispersion. Here expect the variations in pipe 
curvature and torsion to distort any normal profile. Hence expect such vari- 
ations to have a large effect on skewness. 
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Figure 7: coefficients of the dispersion coefficient D for water at 15°C given 
by (33). For higher Reynolds number the highest powers in 1Z dominate the 
coefficients. 



The third order modification to the Taylor model of dispersion is 

d 2 C* 



dt ds ds V ds J ds 



E- 



where the skewness coefficient 



E 



V 3 



2880 



-p 2 ro>6 / 3241338875 4 , 1104359125 2\ _ 4085615 n>4 2 

' L 1 107602919424 35867639808 / 20901888 

, ny2 ( 68855 2 , 985 \ _ 13] 

^ V 96768 12096/ 32 J 

I K J [o>8 ( 5943982203125 6 _ 1183708683125 5 _ 12012557065375 4 

L V 9038645231616 2259661307904 ° 9038645231616 
_ 4304114889625 3\ 
5931610933248° I 



(34) 



(35a) 



(35b) 



■ 2 r 2 + «' 



./2 



432" / ~ 12 J 
[(Y>8 / 254564364353125 6 _ 9130311425570375 4 
L V 1434131710083072 34419161041993728 



(35c) 
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_ 545957180588375 2\ 
29041167129182208 ) 
, <v>6 ( 9773515705025 4 , 241252377733475 2\ 
t -A ^ 2259661307904 759246199455744 ) 

_|_ m 4 / 1227346555 4 _ 24438946685 2 _ 159535 \ 
V 114960384 16554295296 2128896/ 

, o>2 / 1325215 2 _ 118243 \ , _801_1 f^rH 
' ^ \ 258048 580608/ ' 2560 J l OJU i 

-I- V(kk'Y I^R 8 f 26098554271144375 6 , 10015284615625 5 
-|- / ^fi, ) \_J\ y 206514966251962368 ^ 204875958583296 

, 83884205830882375 4 , 68444215116292625 3 , 44871973001125 2\ 
' 206514966251962368 ' 464658674066915328 ' 691456360218624 / 
, o>6 ( 71359651375 5 _ 516258977875 4 _ 596762793925 3 
V 6025763487744 111588212736 2875932573696 
_ 2958969842375 2 _ 62332574575 \ 
23007460589568 1977203644416 / 

, (y>4 ( 48287105 4 _ 480287515 3 , 4233379915 2 , 89672375 , 159535 \ 
' V 3784704 172440576 2759049216 459841536 ' 2128896/ 

, (v>2 / 23023685 2 , 2964127 , 501107 \ _ 22571 , n( 4 t3n /or \ 
^ JX \ 4644864 U T 4644864 U ^ 2322432 J 5760 J ^ ^ ' / ' 

The skewness coefficient for flow in a straight pipe, — "P 3 /2880 from (35a), 
is well known [5]. One outstanding puzzle in the field of dispersion is that 
in rivers one observes contaminant concentrations with long tails upstream 
(not downstream) [2, e.g.]. But theoretical models predict either only a weak 
enhancement of upstream tails or more confoundedly, as the above nega- 
tive skewness coefficient implies for straight pipe flow, a weak downstream 
tail. However, as we now see, curvature effects, presumably induced by the 
secondary flows, significantly change the skewness coefficient thereby enhanc- 
ing the upstream tail of a contaminant release. With the caveat that this 
derivation is for pipes, not rivers, this is a qualitative improvement in the 
theoretical model compared with observations. 



For large Schmidt number cr (typical for the dispersion of material in 
liquids), and similar to the dispersion coefficient D, there are two distin- 
guished limits of the above expression for the skewness coefficient. 

• Firstly, in terms of the magnitude S of the slow axial variations, the 
appropriate scaling is that the fluid flow is slow, 1Z ~ 5, the Schmidt 
number large enough, a ~ 1/S 2 , so that the Peclet number is also large, 
V ~ 1/5, then the skewness coefficient is large, E ~ 1/5 3 . Recalling 
the parameter a = VSP /100 = 7£ 2 cr/100 and evaluating fractions, the 
leading order terms in the skewness coefficient (35) are 

V 3 

E « \1 - k 2 (.8675 a 2 + 20.49) 

2880 L v ' 

+ Vkk' (.1894 a 2 - .2610 a + 11.32) 

+ V 2 (k 2 t 2 + k' 2 ) (-.05112 a 2 + 3.075) 

+ V 2 (kkJ (.03640 a 2 - .003411 a - 3.674)] (36) 
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• Secondly, for a typical Schmidt number a bigger than 10 3 or so, and 
for any flow with Reynolds number ft bigger than about 2, then the 
parameter a will be bigger than about 40 and the above skewness 
coefficient (36) is dominated by the quadratic powers in a: 



In this regime, even small curvature, through the k 2 term, will cause 
the skewness coefficient to become positive, possibly large, and so lead 
to concentration tails upstream (qualitatively as observed in rivers). 
Torsion in the pipe leads to the same upstream tails. 

Recognise another upstream memory effect. The subexpression 0.8675 k 2 — 
0.1894 Vkk' appearing in the first line of (37) is equivalent to simply 0.8675 k 2 
evaluated at a distance £ = 0.1092 V upstream from any particular location. 
This upstream memory is approximately half that of the dispersion coeffi- 
cient. 

For the dispersion of heat in water with a Prandtl number of cr = 8.1 
at 15°C the skewness coefficient (35) reduces to 

E = -0.1845 ft 3 

+ k 2 (-3.291 ft + 3.788 ft 3 - .01039 ft 5 + .001067ft 7 ) 

+ kk (.5833 + 9.956 TZ 2 - 16.43 TZ 4 + .08625 TZ 6 - .002101 TZ 8 ) 

+ (k 2 t 2 + k 2 ) (2.534ft + 27.28ft 3 - 37.30ft 5 + .1509ft 7 + .003968ft 9 ) 

+ (««')' (-3.174 ft - 25.91 ft 3 + 43.37ft 5 - .1589 ft 7 - .002604 ft 9 ) .(38) 

For Reynolds number ft > 11 the highest powers in ft dominate. Observe 
that for these Reynolds numbers both curvature and torsion may easily re- 
verse the sign of of the skewness paramater E through the combination 



Again this effect promotes upstream tails in the dispersion. Although the 
terms in kk' and (/?«/)' may keep E negative, we prefer to interpret these as 
representing upstream memory. 



E 




(37) 



+ft 7 .001067 k 2 + .003968 1Z 2 (k 2 t 2 + k' 2 ) 
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For the dispersion of heat in air with a Prandtl number of a = 0.71 
at 15°C and recalling that m = 7?./10, the dispersion coefficient (35) reduces 
to 

E = -0.1242 m 3 

+ k 2 (-2.884 m + 3. 124 m 3 - 0.6995 m 5 + 0.1645 m 7 ) 

+ kk! (0.5833 + 6.2171H 2 - 9.592 m 4 + 3.537m 6 - 0.7762 m 8 ) 

+ (k 2 t 2 + k' 2 ) (2.221 m + 16.93 m 3 - 25.07 m 5 + 8.940 m 7 - 0.3844 m 9 ) 

+ (kk')' (-2.782 m- 12.99 m 3 + 22.94 m 5 - 9.478 m 7 + 1.287 m 9 ) . (39) 

For Reynolds number 1Z > 40 the highest powers in 1Z dominate. Observe 
that for such larger Reynolds number the sign of the skewness coefficient 
changes sign sensitively depending upon the torsion r, curvature k, and its 
gradients. 



3.4 Higher order curvature affects the dispersion 

Computer algebra also straightforwardly determines even higher order cor- 
rections to the dispersion coefficient. These terms may be used, for example, 
to give estimates of the errors in the earlier approximations. However, the 
algebraic expressions quickly become extremely complicated. We just extend 
the analysis to the next order in curvature, but no higher order in gradients, 
to obtain the following correction to the dispersion coefficient (28): 



D 



|k 4 



605<) 156407 



48 



3094629863915520000 
11319036743801 \ 



a 4 + 



148720297230839 
464658671066915328000000 



a 



14297189971289702400000 
21839753491553 



+n 6 ( 
+n 4 (- 



1265 110332 1262 100000 



a 2 + 



1800408289399 



24648813997 ^2 , 



2711593569484800000/ 
5096950451 \ 



64377815040000 
4685593 77 2 , 13829 ] 
' V ' 9216 J 



348364800 



21459271680000 , 
6 



(40) 



or approximately 



D « + W 



7>V 



+ 



18 



n 
To 



(0.22a 4 + 0.032a 2 - 0.079) 



(^) 6 (l-7- 2 + 0.66) + (^) 4 (-3.8a 2 + 2.4) 
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Figure 8: comparison of the Pade approximant (42) (solid) with experimen- 
tal estimates (circles) collated by Johnson [14], Fig. 9, and the predictions 
(dashed line) of their spectral method based on Poiseuille flow. 



- 1.3 



1.5 



0{k\5). 



(41) 



I also computed the dispersion coefficient to the next correction, with er- 
rors 0(k 8 ,5). Then recalling the Dean number D = the dominant 
terms for the coefficient of dispersion of material in liquids, large Schmidt 
number cr, are: 



D 



^2 

48 



0.3754 



0(5> 16 <7* 



gV 
400 



0.2249 



gv 

400 



-0.1388 



gv 

400 



This expression is valid for small enough D 2 a . Using the additional informa- 
tion [14, §4.3] that the limit at large D 2 a is approximately 0.20 , I construct 
the following Pade approximant in terms of a = D 2 a/ 400 



D „ V2 x 1 + 0.3068 a 2 - 
~ 48 X 1 + 0.6822 a 2 



0.007811a 4 



0.03905 a 4 



(42) 



Tony Roberts, February 8, 2008 



4 Conclusion 



26 



See in Figure 8 that this expression for the dispersion coefficient matches 
reasonably well with experiments over the whole range of D 2 a . 

4 Conclusion 

Computer algebra handles the considerable details of deriving the compli- 
cated expressions describing dispersion in generally curving pipes, §A. Fix- 
ing one of the fluid flux or the mean pressure gradient affects the dispersion, 
§1, and throughout I present results for the appropriate case of fixed fluid 
flux. The Pade approximation (42) for the dispersion in a constant curvature 
pipe is reasonably accurate over the entire range of Dean numbers. When 
the pipe's curvature and torsion vary, much of the effects of variations upon 
the dispersion may be recast as an upstream memory. Overall, torsion r in 
the pipe seems to have little effect on the dynamics except for a sensitivity, in 
combination with the curvature k, of the skewness, §3.3. The skewness coef- 
ficient is very sensitive to curvature, hence is easily made positive, and which 
may thus explain the observations of long upstream tails in the dispersion of 
material in rivers. 

A Computer algebra derivation 

Computer algebra is a very powerful means to derive asymptotic expansions. 
In order for other people to reproduce and verify the results recorded herein, 
I list here the core of the program used to derive the asymptotic expansions; 
obtain the full program by request. 

The computer algebra program was written in reduce 5 to calculate the 
asymptotic expansions of the centre manifold models described in this article. 

There is a lot of detail to the computer algebra program. However, the 
key to the correctness of the results is the coding of the governing equations 
which forms the key part of the core printed here. The algorithm iteratively 
drives to zero the residuals of these equations, see Roberts [26] for a generic 
description of the algorithm. Thus the details about how the residuals are 
reduced are not vital, only that they are correctly computed and ultimately 
zero. 

5 At the time of writing, information about reduce was available from Anthony 
C. Hcarn, RAND. Santa Monica, CA 90407-2138, USA. mailto:reduce@rand.org 
There were demonstration versions of reduce freely available at 
f tp : //f tp . zib . de/pub/ reduce/demo or f tp : //ftp . maths . bath . ac . uk/pub/ algebra . 
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1 comment Find the flow in an arbitrarily curving pipe, 

2 Simultaneously determine the shear dispersion in such a flow: 

3 eps=magnitude of curvature terms, 

4 del=magnitude of axial derivatives and of torsion. 

5 ; 

6 depend kap , s ; % curvature 

7 depend tau,s; % torsion 

8 % local coordinate system is (s,r,th), tp=th+phi(s) 

9 depend tp,s,th; 

10 let { df (tp,s)=>-tau, df(tp,th)=>l >; 

11 hs : =l-eps*kap*r*cos (tp) ; 

12 hr:=l; 

13 ht:=r; 

14 % trigonometry rules OK 

15 let { sin(~a)*cos(~b) => (sin(a+b)+sin(a-b) )/2 



16 , cos(~a)*cos(~b) => (cos(a-b)+cos(a+b))/2 

17 , sin(~a)*sin(~b) => (cos (a-b) -cos (a+b) ) /2 

18 , cos(~a)~2 => (l+cos(2*a))/2 

19 , sin(~a)~2 => (l-cos(2*a))/2 

20 }; 



21 % mean over a cross section (mult by r to use) 

22 depend r,rt; depend tp,rt; 

23 operator mean; linear mean; 

24 let { mean(r~~m*cos (~n) ,rt) => 



25 , mean(r~~m*sin(~n) ,rt) => 

26 , mean(r~~m,rt) => 2/(m+l) 

27 , mean(r,rt) => 1 

28 }; 



29 % operators to solve for updates 

30 ... 

31 % initial approximations 

32 u:=2*(l-r~2) +eps*3/2*kap* (r-r~3) *cos (tp) ; 

33 v:=0; 

34 w:=0; 

35 % pressure = ps + p 

36 % = (local mean gradient) + (zero mean fluctuation) 

37 ps:=-8; p:=0; 

38 % concentration of tracer, mean c. 

39 depend c,s,t; 

40 let df(c,t)=>g; 

41 cc:=c; 

42 g:=0; 

43 pe:=re*sc; % Peclet number = Reynolds * Schmidt 

44 rh:=l; °L approx reciprocal of axial scale factor hs 
45 

46 7, iterate until residuals are negligible 

47 let { eps~3=>0, del~5=>0 }; °/« truncate the asymptotics 

48 repeat begin 

49 '/„ reciprocal of scale factor 



Tony Roberts, February 8, 2008 



A Computer algebra derivation 



28 



50 eqr : =hs*rh-l ; 

51 rh:=rh-eqr; 

52 7o vorticity 

53 oms:=(df (r*w,r)-df (v,th))/r; 

54 omr : = (df (hs*u , th) -del*df (r*w , s) ) *rh/r ; 

55 omt : = (del*df (v , s) -df (hs*u , r ) ) *rh ; 

56 °/ Navier-Stokes equation 

57 nss:=rh*(ps+del*df (p,s)) +(df (r*omt,r)-df (omr,th))/r 

58 +re*( del*u*df (u, s) *rh+v*df (u,r)+w*df (u,th)/r 

59 +u*df (hs,r)*v*rh+u*df (hs,th)*w*rh/r ); 

60 nsr:=df(p,r) +(df (hs*oms ,th) -del*df (r*omt , s) ) *rh/r 

61 +re*( del*u*df (v, s) *rh+v*df (v,r)+w*df (v,th)/r 

62 -w~2/r-df (hs,r)*u~2*rh ) ; 

63 nst:=df (p,th)/r +(del*df (omr , s) -df (hs*oms ,r) ) *rh 

64 +re*( del*u*df (w, s) *rh+v*df (w,r)+w*df (w, th)/r 

65 -u~2*df (hs,th)*rh/r+v*w/r ) ; 

66 °/ continuity equation 

67 cty : = (del*df (r*u , s) +df (r*hs*v , r ) +df (hs*w , th) ) *rh/r ; 

68 ... 

69 equation for tracer evolution 

70 ceq:= df(cc,t) +pe*( del*u*df (cc , s) *rh+v*df (cc ,r)+w*df (cc , th) /r ) 

71 -(del"2*df (rh*r*df (cc,s) ,s)+df (r*hs*df (cc,r) ,r) 

72 +df (hs/r*df (cc.th) ,th))*rh/r; 

73 cmean:=mean(r*hs*cc,rt)-c; 
74 

75 end until (eqr=0)and(nss=0)and(nsr=0)and(nst=0) 

76 and(cty=0)and(ceq=0)and(cmean=0) ; 



77 

78 % check subsiduary conditions 

79 uwall:=sub(r=l,u) ; 

80 vwall:=sub(r=l,v) ; 

81 wwall:=sub(r=l,w) ; 

82 umean:=mean(r*u,rt) ; 

83 pmean:=mean(r*p,rt) ; 

84 cwall:=sub(r=l,df (cc,r)) ; 

85 cf lux : =mean(r* (u*pe*cc-del*rh*df (cc , s) ) ,rt) ; 

86 end; 

Observe that the pressure is decomposed into a mean gradient and a 
cross-pipe fluctuating component. There is code to adjust the mean pressure 
gradient to ensure a constant mean fluid flux. To adapt to the traditional 
fixed pressure gradient in a helical or toroidal pipe, one just needs to omit 
the modifications. However, the results are then inappropriate to a pipe with 
varying curvature or torsion. 



Tony Roberts, February 8, 2008 



References 



29 



References 

[1] G. K. Batchelor. An Introduction to Fluid Dynamics. CUP, 1979. 3, 8 

[2] T. Beer and P. C. Toung. Longitudinal dispersion in natural streams. 
J. Environmental Eng, 109:1049-1067, 1983. 22 

[3] S. A. Berger, L. Talbot, and L.-S. Yao. Flow in curved pipes. Ann. Rev. 
Fluid Mech., 15:461-512, 1983. 3, 4, 5, 12, 14 

[4] P. C. Chatwin. The approach to normality of the concentration distri- 
bution of a solute in a solvent flowing along a straight pipe. J. Fluid 
Mech, 43:321-352, 1970. 20 

[5] P. C. Chatwin and C. M. Allen. Mathematical models of dispersion in 
rivers and estuaries. Annu. Rev. Fluid Mech., 17:119-149, 1985. 22 

[6] P. Daskopoulos and A. M. Lenhoff. Flow in curved ducts: bifurcation 
structure for stationary ducts. J. Fluid Mech., 203:125-148, 1989. 3 

[7] W. R. Dean. Note on the motion of fluid in a curved pipe. Phil. Mag., 
4:208-223, 1927. 6 

[8] W. R. Dean. The stream-line motion of fluid in a curved pipe. Phil. 
Mag., 5:673-695, 1928. 6 

[9] M. E. Erdogan and P. C. Chatwin. The effects of curvature and buoyancy 
on the laminar dispersion of solute in a horizontal tube. J. Fluid Mech., 
29:465-484, 1967. 6 

[10] D. Fargie and B. W. Martin. Developing flow in a pipe of circular cross 
section. Proc. Roy. Soc. Lond. A, 321:461-476, 1971. 14 

[11] D. Gammack and P. E. Hydon. Flow in pipes with non-uniform curva- 
ture and torsion. J. Fluid Mech., 433:357-382, 2001. 3, 6, 9, 13 

[12] M. Germano. On the effect of torsion on a helical pipe flow. J. Fluid 
Mech., 125:1-8, 1982. 3, 7 

[13] Sandip Ghosal. Lubrication theory for electro-osmotic flow in a microflu- 
idic channel of slowly varying cross-section and wall charge. J. Fluid 
Mech., 459:103-128, 2002. 4 

[14] M. Johnson and R. D. Kamm. Numerical studies of steady flow dis- 
persion at low Dean number in a gently curving tube. J. Fluid Mech., 
172:329-345, 1986. 3, 5, 12, 25 



Tony Roberts, February 8, 2008 



References 



30 



[15] S. W. Jones and W. R. Young. Shear dispersion and anomalous diffusion 
by chaotic advection. J. Fluid Meek, 280:149-172, 1994. 3 

[16] H. C. Kao. Torsion effects on fully developed flow in a helical pipe. 

J. Fluid Meek, 184:335-356, 1987. 7, 12 

[17] E. Kreyszig. Advanced engineering mathematics. Wiley, 8th edition, 
1999. 4 

[18] J. Larrain and C. F. Bonilla. Theoretical analysis of pressure drop in 
the laminar flow of fluid in a coiled pipe. Trans. Soc. Rheoi, 14:135-147, 
1970. 10 

[19] S. Liu and J. Masliyah. Axially invariant laminar flow in helical pipes 
with finite pitch. J. Fluid Meek, 251:315-353, 1993. 3, 7 

[20] D. J. McConalogue. The effects of secondary flow on the laminar disper- 
sion of an injected substance in a curved tube. Proc. Roy. Soc. Lond. A, 
315:99-113, 1970. 3 

[21] G. N. Mercer and A. J. Roberts. A complete model of shear dispersion 
in pipes. Jap. J. Indust. Appl. Math., 11:499-521, 1994. 6 

[22] S. Murata, Y. Miyake, and T. Inaba. Laminar flow in a curved pipe 
with varying curvature. J. Fluid Mech., 73:735-752, 1976. 7, 12, 14 

[23] R. J. Nunge, T.-S. Lin, and W. N. Gill. Laminar dispersion in curved 
tubes and channels. J. Fluid Mech., 51:363-383, 1972. 3, 5, 6 

[24] T. J. Pedley. Fluid mechanics of large blood vessels. C.U.P., 1980. 3, 
7, 12 

[25] A. J. Roberts. The application of centre manifold theory to the evolution 
of systems which vary slowly in space. J. Austral. Math. Soc. B, 29:480- 
500, 1988. 4 

[26] A. J. Roberts. Low- dimensional modelling of dynamics via computer 
algebra. Computer Phys. Comm., 100:215-230, 1997. 26 

[27] A. J. Roberts. Low- dimensional modelling of dynamical systems applied 
to some dissipative fluid mechanics. In Rowena Ball and Nail Akhme- 
diev, editors, Nonlinear dynamics from lasers to butterflies, volume 1 
of Lecture Notes in Complex Systems, chapter 7, pages 257-313. World 
Scientific, 2003. 4 



Tony Roberts, February 8, 2008 



References 



31 



[28] D. M. Ruthven. The residence time distribution for ideal laminar flow 
in helical tubes. Chem Eng Sci, 26:1113-1121, 1971. 3 

[29] R. Smith. Longitudinal dispersion coefficients for varying channels. 
J. Fluid Meek, 130:299-314, 1983. 6, 18 

[30] G. I. Taylor. Dispersion of soluble matter in solvent flowing slowly 
through a tube. Proc. Roy. Soc. Lond. A, 219:186-203, 1953. 4 

[31] G. I. Taylor. Conditions under which dispersion of a solute in a stream 
of solvent can be used to measure molecular diffusion. Proc. Roy. Soc. 
Lond. A, 225:473-477, 1954. 4 

[32] H. C. Topakoglu. Steady laminar flows of an incompressible viscous fluid 
in curved pipes. J. Math. Mech., 16:1321-1338, 1967. 6 

[33] R. N. Trivedi and K. Vasudeva. Axial dispersion in laminar flow in 
helical coils. Chem Engrg. Sci., 30:317-325, 1975. 3, 4, 5 

[34] E. R. Tuttle. Laminar flow in twisted pipes. J. Fluid Mech., 219:545-570, 
1990. 3, 4, 8, 9, 12 

[35] S. D. Watt and A. J. Roberts. The construction of zonal models of 
dispersion in channels via matching centre manifolds. J. Austral. Maths. 
Soc. B, 38:101-125, 1996. 

[36] K. Yamamoto, S. Yanase, and T. Yoshida. Torsion effect on the flow in 
a helical pipe. Fluid Dyn. Res., 14:259-273, 1994. 6, 7 

[37] K. Yamamoto, S. Yanase, and T. Yoshida. Erratum: Torsion effect on 
the flow in a helical pipe. Fluid Dyn. Res., 24:309-311, 1999. 6, 7 

[38] L. Zabielski and A. J. Mestel. Steady flow in a helically symmetric pipe. 
J. Fluid Mech., 370:297-320, 1998. 3, 9 



Tony Roberts, February 8, 2008 



